11-4 ODE 仵蛌i峈k

在前一節中，我們已知道基本 ODE 檔案的寫法，如 vdp1.m、vdp2.m 及 lorenzOde.m 等。在本節中，將更進一步介紹 ODE 檔案的進階用法，以使 ODE 指令能夠迅速且準確地算出積分結果。

首先，我們可將 tspan（積分時間範圍）、y0（起始值）及 options（ODE參數）置於 ODE 檔案中，這些變數必須能由 ODE 檔案傳回，其格式為：

[tspan, y0, options] = odeFile([], [], 'init')

其中我們假設 odeFile 即是我們的 ODE 檔案。若 odeFile 滿足上述要求，則我們可以直接呼叫 ODE 指令如下：

[t, y] = solver('odeFile')

其中 solver 為前述的任一個 ODE 指令，它可由 odeFile 直接得到 tspan、y0 及 options 等積分所需的資訊。

以前述的 van der Pol 為例，若要能夠傳回 tspan、y0 及 options，vdp1.m 須改寫如下（vdp3.m）：

11-常微分方程式/vdp3.m

function [output1, output2, output3] = vdp3(t, y, flag) if strcmp(flag, '') mu = 1; output1 = [y(2); mu*(1-y(1)^2)*y(2)-y(1)]; % dy/dt elseif strcmp(flag, 'init'), output1 = [0; 25]; % Time span output2 = [3; 3]; % Initial conditions output3 = odeset; % ODE options end

此時我們可以直接呼叫 ODE 指令如下：

Example 1: 11-常微分方程式/odeAdvanced01.m

ode45('vdp3')

此外，van der Pol 的微分方程式有一個參數 $\mu$，若希望從外面傳入此參數的值，我們可將 vdp3.m 改寫成vdp4.m，其內容如下：

11-常微分方程式/vdp4.m

function [output1, output2, output3] = vdp4(t, y, flag, mu) if nargin < 4 | isempty(mu), mu = 1; end if strcmp(flag, '') output1 = [y(2); mu*(1-y(1)^2)*y(2)-y(1)]; % dy/dt elseif strcmp(flag, 'init'), output1 = [0; 25]; % Time span output2 = [3; 3]; % Initial conditions output3 = odeset; % ODE options end

此時 $\mu$ 就變成一個選擇性（Optional）的參數，其預設值為 1。在下例中，我們將 $\mu$ 的值從 MATLAB 傳入，並畫出不同 $\mu$ 值下的 van der Pol 方程式的狀態變數：

Example 2: 11-常微分方程式/odeAdvanced02.m

subplot(2,1,1); ode45('vdp4', [], [], [], 1); % mu=1 subplot(2,1,2); ode45('vdp4', [], [], [], 3); % mu=3

在上圖中，$\mu$ 的值分別是 1 及 3。指令列中用到了許多空矩陣，這些空矩陣代表「取用預設值」，因此 ode45 會直接從 vdp4.m 取用時間區間及變數起始值。事實上，您也可以傳入二個或更多的參數，MATLAB 及 ODE 指令對於可傳入的參數個數並無設限。

此外，為解決其它較複雜的 ODEs 及 DAEs（Differential Algebra Equations），ODE 檔案亦可在不同的旗號（Flag）下傳回不同的資訊，以下是一個完整的列表：

整理：由 odeset 產生的 ODE 選項

旗號	傳回值
（空字串）	dy(=F(t,y))
init	tspan, yo 及 options
jacobian	Jacobian 矩陣 J(t,y) = $\partial F/\partial Y$
jpattern	Jacobian sparsity pattern 之矩陣
mass	解 M(t,y)y' = F(t,y) 所須的質量矩陣 M
events	定義事件發生點的各種資訊

後四者因為不常用到，在此不再贅述。有興趣的讀者可參考 MATLAB 有關 ODE 的手冊。

---

MATLAB程式設計：進階篇